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We study the melting transition of the low-temperature vortex solid in strongly anisotropic layered 
superconductors with a concentration of random columnar pinning centers small enough so that the 
areal density of the pins is much less than that of the vortex lines. Both the external magnetic field 
and the columnar pins are assumed to be oriented perpendicular to the layers Our method, involving 
numerical minimization of a model free energy functional, yields not only the free energy values at 
the local minima of the functional but also the detailed density distribution of the system at each 
minimum: this allows us to study in detail the structure of the different phases. We find that at 
these pin concentrations and low temperatures, the thermodynamically stable state is a topologically 
ordered Bragg glass. This nearly crystalline state melts into an interstitial liquid (a liquid in which a 
small fraction of vortex lines remain localized at the pinning centers) in two steps, so that the Bragg 
glass and the liquid are separated by a narrow phase that we identify from analysis of its density 
structure as a polycrystalline Bose glass. Both the Bragg glass to Bose glass and the Bose glass to 
interstitial liquid transitions are first-order. We also find that a local melting temperature defined 
using a criterion based on the degree of localization of the vortex lines exhibits spatial variations 
similar to those observed in recent experiments. 

PACS numbers: 74.25.Qt, 74.72.Hs, 74.25.Ha, 74.78.Bz 



I. INTRODUCTION 

The mixed phase of type-II superconductors with ran- 
dom pinning is generally regarded to be an archetypal 
test system for the study of the effects of quenched disor- 
der on the structure and melting of solids. In this phase, 
magnetic flux penetrates the sample as quantized vortex 
lines vifhich form a special physical system known as "vor- 
tex matter" . The fascinating equilibrium and dynamical 
properties of vortex matter in the mixed phase of high- 
temperature superconductors (HTSCs) have prompted 
considerable experimental and theoretical attention for 
more than a decade (see Ref. for an early review). The 
mixed phase of HTSCs is strongly affected by the pres- 
ence of pinning centers, either intrinsic or artificially gen- 
erated. Understanding the effects of pinning in these sys- 
tems is basic for practical applications because pinning 
strongly influences the value of the critical currenli^. 

The vortex lines in a pure type-II superconductor form 
a triangular lattice (the Abrikosov lattice) at low tem- 
peratures. Because of enhanced thermal fluctuations in 
highly anisotropic, layered HTSCs, the Abrikosov lattice 
in very pure samples, in a magnetic field perpendicular 
to the layers (the field will be assumed to be in this direc- 
tion throughout our discussion), undergoes a first-order 
melting transition^ into a resistive vortex liquid (VL) as 
the temperature T increases. Random pinning destroys 
the long-range translational order of the Abrikosov lattice 
and leads to the occurrence of a variety of glassy phases 
at low T. It is now generally accepted that in systems 



with random point pinning, a topologically ordered low- 
temperature phase with quasi-long-range translational 
order (denoted as the "Bragg glass" (BrG) phase) occurs 
at low fields if the pinning disorder is sufficiently weak. 
This has been established theoreticalljiSi^ as well as ex- 
perimentally (see e.g. Ref.'^. The possibility of an amor- 
phous vortex glass (VG) phase, with nonlinear voltage- 
current characteristics and vanishing resistance in the 
zero-current limit, in systems with strong pinning (or at 
high magnetic fields where the effects of pinning disorder 
are enhanced) was suggested by Fishe»2, (see also Ref. 8) . 
However, in spite of extensive investigations, the exis- 
tence of a true VG phase (i.e., an amorphous glassy phase 
thermodynamically distinct from the high-temperature 
VL) in systems with uncorrelated point pinning remains 
very controversial: different calculations2ii£ lead to dif- 
ferent conclusions and the experimental situationiii^ is 
similarly contradictory. A variety of "glassy" behavior 
has been reported in different experiments (see Ref. 
and references therein) near the first-order melting tran- 
sition of the BrG phase of both conventional supercon- 
ductors and HTSCs with point pinning. It has been 
suggested^^'^** that these observations may be understood 
if it is assumed that the melting of the BrG phase occurs 
in two steps: the BrG first transforms into a "multido- 
main" glassy phase that melts into the usual VL at a 
slightly higher T. 

Columnar pinning defects can be produced for example 
as damage tracks arising from heavy- ion bombardment. 
The technological importance of these extended defects 



2 



oriented parallel to the direction of the external magnetic 
field has been long recognizediSii^: they are highly effec- 
tive in increasing the critical current by localizing vortex 
lines along their length. Heavy-ion irradiation produces a 
random array of parallel columnar defects each of which 
can trap one or more vortices at low temperatures. The 
effects of such random arrays of extended defects, ori- 
ented perpendicular to the superconducting layers, on the 
properties of the mixed phase of HTSCs have been ex- 
tensively studied experimentalljiiiiSiiS. The same ques- 
tion has also been examined theoretically22i2i*22i2^ and 
numericallji2ii2Si2S42i. When the columnar pinning is 
strong, and their concentration exceeds that of the vor- 
tex lines, a so-called "strong" Bose Glass (BoG) phase 
with nearly all the vortex lines localized at the pinning 
centers occurs^*''^'* at low T. This phase is strongly dis- 
ordered with very short-range translational and bond- 
orientational correlations. The behavior in this regime 
is fairly well-understood in terms of a mapping^*^ of the 
thermodynamics of the system to the quantum mechani- 
cal properties of a two-dimensional system of interacting 
bosons in a random potential. 

Much less is known about the behavior in the dilute 
pin limit where the concentration of pins is much smaller 
than that of vortex lines. In such systems, one expects^^ 
a "weak" BoG phase at low temperatures, with a small 
fraction of vortex lines localized strongly at the pinning 
centers and the remaining ones localized relatively weakly 
in the interstitial region between pinning centers. As the 
temperature is increased, this phase shoul d? melt into 
an interstitial liquid (IL). In the IL phase, some of the 
vortices remain trapped at the pinning centers, while the 
other, interstitial ones, form a liquid. The pinned vor- 
tices are expected2i*2^ to delocalize, thereby forming the 
usual VL, at a crossover occuring at a higher tempera- 
ture. The melting transition of the "weak" BoG phase 
into the IL is predictedSiiS^ to be first-order for small pin 
concentrations, whereas the "strong" BoG to VL transi- 
tion occurring for large pin concentrations is knownSSiSS 
to be a continuous one. 

A first-order melting transition of the "weak" BoG 
phase has been observed^ in recent experiments on sam- 
ples of Bi2Sr2CaCu208+x (BSCCO) with a small con- 
centration of columnar pins. The BoG phase in these 
systems is foundiSiiS to have a polycrystalline struc- 
ture with ordered vortex crystallites of different orien- 
tations embedded in the interstitial region between vor- 
tices pinned at the columnar defects. If the pin con- 
centration is sufficiently small, the melting of this BoG 
phase upon increasing the temperature occurs into an IL 
phase with a fraction of vortex lines remaining pinned at 
the defects. These pinned vortices delocalize at a higher 
temperaturSiiS. The temperature at which the first-order 
BoG to IL transition occurs isi^ii^ very close to the melt- 
ing temperature of the same system without the colum- 
nar pins, if the pin concentration is small. As the pin 
concentration is increased, the melting temperature in- 
creases and the transition eventually becomes continu- 



ous. The difference between the melting temperature 
and the temperature at which the pinned vortices in 
the IL phase delocalizes decreases and eventually goes 
to zero as the pin concentration increases. Another in- 
teresting feature, found experimentally for both point^ 
and columnar— pinning, is that the melting of the low 
T "solid" phase is "inhomogeneous" in that it occurs 
locally over a range of temperatures: the local transi- 
tion temperature, which can be measured from the lo- 
cal magnetization, is different in different regions of the 
sample. This inhomogeneity of the local melting tem- 
perature is believed^^ to be closely related to the local 
arrangement of the specific pinning centers (the so-called 
"pinning landscape") in each sample studied. 

The polycrystalline nature of the BoG phase for 
small pin concentrations has also been observed in 
simulations^^. Recent numerical work^i indicates that 
a BrG phase, with topological order, may also exist in 
such systems provided that the pin concentration is suf- 
ficiently small. In these simulations, the BrG phase is 
found to melt into the IL phase via a first-order tran- 
sition as the temperature is increased. As the pin con- 
centration is increased beyond a critical value, the BrG 
phase disappears^^ and a low-temperature BoG phase 
with a continuous BoG to IL transition upon increasing 
the temperature is foundSSiSi. 

In this paper, we approach these problems through 
a different numerical method. We consider a layered, 
strongly anisotropic superconductor (such as BSCCO) 
with a dilute random array of columnar defects oriented 
perpendicular to the layers, and with a magnetic field 
applied parallel to the columnar pins. We describe the 
equilibrium properties of this system in terms of vor- 
tex density variables, using a free energy functional of 
the Ramakrishnan-Yussouff^^ (RY) form. We numeri- 
call minimize a spatially discrctized version of this free 
energy functional and obtain the vortex density configu- 
ration at each local minimum. Analysis of these density 
configurations, both in terms of a variety of correlation 
functions and by direct visualization of the arrangement 
of the vortices (vortex positions are identified with those 
of local peaks of the vortex density), allows us to identify 
the nature of the phases corresponding to different local 
minima of the free energy. Comparison of the values of 
the free energy at the minima as the temperature varies 
yields the transition temperatures. A similar procedure 
has been successfully employed for the case of a regular 
array22i^ of columnar pinning centers. 

We find in our study that that a Bragg glass phase 
exists at low temperatures for samples with a small con- 
centration of columnar pins. We will present evidence for 
this phase from our analysis of the vortex density and as- 
sociated correlation functions. Upon warming, this phase 
melts into what we show to be an insterstitial liquid (IL) 
phase, but the melting is shown^ to occur in two steps: 
the BrG and IL phases are separated, over a narrow tem- 
perature range, by an intermediate phase which exhibits 
a multi-domain polycrystalline structure. We show that 
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this intermediate phase should be identified as a Bose 
glass (BoG). The temperature of the upper (BoG to IL) 
transition is approximately independent of the columnar 
pin concentration c, at the low c values studied. Both the 
BrG to BoG and the BoG to IL transitions arc found to 
be first-order. We also find that a local melting temper- 
ature can be suitably defined using a criterion based on 
the degree of localization of the vortices, and that its be- 
havior (in particular, its spatial variation) is quite consis- 
tent with what is seen in experimentsiStS^. The value of 
the local melting temperature is strongly correlated with 
the presence of topological defects (dislocations) in the 
vortex solid which, in turn, is correlated with the local ar- 
rangement of the pinning centers. We also show that the 
transition to the IL phase corresponds to a percolation of 
regions containing liquid-like (delocalized) vortices across 
the sample. 

After this Introduction, we discuss our definitions, 
model and numerical methods in the next Section. Then 
in Section IIIII we present our results. We discuss first 
the free energy minima and their study through the use 
of tools such as correlation functions, local peak-density 
plots, and Voronoi construction for the lattice formed by 
local density peaks. These tools allow the identification 
of the phases at each free energy minimum, as we shall 
show. Then we derive the phase diagram and show that 
indeed the BrG and IL phases are separated by a thin 
sliver of BoG. Finally, our results for the nature of the 
local melting are displayed and discussed. A brief conclu- 
sions Section recapitulates the main points of the paper 
and discusses them in the context of existing results. 

II. METHODS 

The general procedures that we use are quite similar 
to those employed in previous work on a regular array 
of columnar pin&22j2i. We will therefore give here only a 
brief summary, emphasizing the details that are different 
in the random case considered here. 

The system we study is a layered superconductor in the 
extreme anisotropic limit, that is, with vanishing Joseph- 
son coupling between layers, which are then coupled via 
the electromagnetic interaction only. This limit is appro- 
priate for the Bi- and Tl- based HTSG compounds in a 
large region of the magnetic field {H)-T plane. In this 
work, we will use material parameter values appropriate 
to BSGCO. With these assumptions, one can write the 
energy of a system of "pancake" vortices residing in the 
layers as a sum of anisotropic two-body interactions. For 
straight columnar pins normal to the layers, the pinning 
potential is the same in every layer. It is then possible to 
write the free energy as a functional of the time averaged 
areal vortex density p{r): 

F[p]-Fo=FRY[p]+Fp[p]. (2.1) 

where Fq is the free energy corresponding to a uniform 
vortex liquid of density po = i?/$o (B is the magnetic in- 



duction and $0 the superconducting flux quantum). The 
first term in the right-hand side of Eq. 1)2. l|l is the free en- 
ergy functional in the absence of pinning. As explained 
above we use for this free energy the BY— functional, 
which is known^fli^iiSi^ to give a quantitatively accu- 
rate description of the vortex-lattice melting transition 
in the absence of pinning. It is of the form: 

PFhy[p] = 1 d^r [p(r){ln(p(r)) - In(po)} - Mr)] 

_ 1 y d'^r' C{\r ^ r'\)Sp{r)Spir'l2.2) 

Here f3 is the inverse temperature, Sp{r) = p(r) ~ po, and 
C(r) is the usual direct pair correlation function^S, which 
may be written as a sum over layers: C(r) = C'(n, r), 
with C(n, r) (where n is the layer separation and r the 
in-layer distance) being the corresponding direct pair cor- 
relation function of a layered liquid of pancake vortices. 
The direct correlation function, which is needed as in- 
put in our free energy, can be accurately calculated in 
a number of ways. We will use here the results of the 
hypernetted chain calculation of Ref. Q In general, the 
results in the limit considered depend on the values of 
the in-plane London penetration length A(T), the inter- 
layer spacing d and a dimensionless coupling parameter 
F given by: 

F = pd^ll%'K^}?{T). (2.3) 

For BSCCO we will take d — ISA, and we will assume 
a standard two- fluid temperature dependence for A(r) 
with A(0) = 1500A and = 85K (at zero field). 

The second term in the right-hand side of Eq. 12.1(1 
represents the pinning and is of the form: 

FM = / d\V^(v)8p{v). (2.4) 

where Vp is the pinning potential which can be written 
as 

Fp(r)=^Fo(|r-R,|), (2.5) 
i 

with the sum extending over the planar positions of the 
random pinning centers. We take the potential Vb cor- 
responding to a single pinning center to be of the usual 
truncated parabolic formi^ 

m{r) = -aF[l - (r/ro)2]e(ro - r) (2.6) 

where tq is the range. The basic length in the prob- 
lem, which we will use as our unit of length unless other- 
wise indicated, is ao, defined by the relation 7ra§po — 1- 
We choose ro to be tq = O.lao and take the dimension- 
less constant a = 0.05, for which value, as previously 
showu)^ each pinning center traps approximately one 
vortex in BSCCO, in the temperature range of interest. 
This range is determined by the following considerations. 
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We will keep the field fixed at B = 0.2T, and vary the 
temperature. The melting temperature of the unpinned 
lattice is then T,"^, ~ 18AK^, and we, therefore, consider 
the neighborhood of this temperature. Then one has, to 
a very good approximation T ~ 2650/r (with the tem- 
perature in Kelvin), while A w A(0). 

To carry out our numerical calculations, we discretize 
the density variables on a computational triangular grid 
of lattice spacing h, containing iV^ sites. We take the 
spacing h to he h — a/16 in our calculations, where 
a = 1.998ao is the equilibrium lattice constant of the 
vortex lattice in the temperature range considered at the 
indicated field. This value, as pointed out in Ref. Isil 
is slightly higher than that of a triangular vortex lat- 
tice of density po, which is (27r/-\/3)^/^ao. We define at 
each site j on this computational lattice a variable pj, 
with Pj = p{rj)v, where v is the area of each computa- 
tional unit cell. Results reported here are for N = 1024, 
which, for the chosen value of h, corresponds to includ- 
ing A^^ = 4096 vortices in the calculation. Preliminary 
and confirmatory results at TV = 512 (1024 vortices) were 
also obtained. A number of pinning sites (Np) are put, 
at random, on some of the computational lattice sites. 
The results presented here, which correspond to the di- 
lute limit, will be primarily for a pin concentration c of 
c = 1/64, that is, a number of pins iVp = 64 at iV = 1024 
(or Np = 16 at TV = 512). Some preliminary results 
for c = 1/32 {Np = 128 at iV = 1024, or Np = 32 
at N = 512) will also be discussed briefly. With ran- 
dom pins, the results depend on the particular random 
pin configuration and averaging over different such con- 
figurations is required. The dependence is however not 
strong: as we shall see, averaging over five to ten config- 
urations is enough to make statistical errors sufficiently 
small for our purposes. 

To perform our studies, we numerically minimize the 
discretized free energy with respect to the ^ discretized 
density variables {pi}. To do so, the interaction term 
in the right side of Eq. (|2.2|) must be repeatedly evalu- 
ated, and since this term is of a convolution form, this 
is most efficiently done in momentum space, through the 
use of efficient Fast Fourier Transform (FFT) routines. 
This avoids having to evaluate this term as a double 
sum which would be computationally much more cum- 
bersome than performing the direct and inverse FFT's. 
In performing the minimization, one must keep in mind 
that the variables {pi} must be nonnegative. This pre- 
cludes the use of many efficient minimization algorithms. 
We use a procedure^'' that ensures that this constraint 
is satisfied. Numerical minimization is performed start- 
ing with an appropriate initial condition for the density 
variables. As explained in Sec. IIIII below, different local 
minima may be found, at the same field and tempera- 
ture, depending on the initial conditions. These minima 
correspond to different phases. The minimization pro- 
cedure yields not only the value of the free energy at 
each minimum but also the detailed vortex density con- 
figuration at the minimum found, i.e. the values of the 



set {pi} at the minimum. It is then straightforward to 
analyze the actual density configuration in several ways. 
It is also possible to evaluate any desired density corre- 
lations. The nature of the phase corresponding to each 
local minimum of the free energy can be inferred from 
such analysis. The (mean-field) phase diagram is then 
obtained from a comparison of the free energies of the 
different minima as the temperature T is varied. The 
results of carrying out this program are discussed below. 



III. RESULTS 
A. Free energy minima 

As explained above, minima of the free energy are 
found by starting the minimization process with appro- 
priate initial conditions. We use three different kinds of 
initial conditions. The first kind is a uniform density 
(Pi = Pov for all i), corresponding to a completely disor- 
dered liquid state. This is typically used to obtain, as we 
shall see, liquid-like states at temperatures near the equi- 
librium melting temperature of the pure vortex lattice, 
which ia^Si^i ~ 18.4K for the value of B considered 
here. The second kind corresponds to a crystalline ini- 
tial state, with values of {pi} as obtained in Ref. |^ for 
the pure vortex lattice. In our system with a spatially 
varying random pinning potential, the pinning energy of 
such a crystalline density configuration depends on the 
choice of the computational lattice sites at which the pe- 
riodic local peaks of the density are located. We take, 
among all possible choices equivalent by symmetry op- 
erations in the computational lattice, the one for which 
the pinning energy is lowest for the specific pinning con- 
figuration under consideration. These initial conditions 
are used to obtain ordered states at lower temperatures. 
Finally, states originally obtained by either of these two 
procedures can be slowly warmed up, or cooled down: in 
this case the initial conditions are simply the final state 
obtained at the previous temperature. As one warms 
up or cools down a state, its nature in general changes: 
liquid configurations may become unstable upon cooling 
and ordered configurations upon warming. 

Whichever initial conditions one uses, a local minimum 
of the free energy is eventually found. At a given tem- 
perature, in general several local minima with different 
density configurations, characterized by the values of the 
set of {pi\ variables, are obtained. The minimum with 
the lowest free energy represents, at a given T, the true 
equilibrium state. It is obviously necessary to develop 
systematic procedures to identify the structure of each of 
these diverse minima. Since we have access to the full set 
of density variables at the computational lattice sites, we 
have at our disposal a variety of methods to achieve that 
goal. The first is to calculate the density correlations, 
e.g. the structure factor S'(k): 

^(k) = |p(k)|VA^. (3.1) 



5 



where p(k) is the discrete Fourier transform of the set 
{pi}. Equivalently, one can consider the Fourier trans- 
form of S'(k), which is the two-point spatial correlation 
function g{r) of the time-averaged local density. We will 
present here results for ^(k), considered as a function 
of the vector k, and for the angularly averaged spatial 
density correlation function, g{r). 

It is also very useful to consider the spatial structure 
formed by the vortices at low temperatures. For this 
purpose we have extracted from {pi} the local peak den- 
sities. We say that the density locally peaks at site i if 
the value of pj at j = i is higher than that at any other 
computational lattice site within a distance a/2 from i, 
where, as defined above, a is the equilibrium spacing of 
the unpinned vortex lattice. As expected, we find that 
at low-temperature solid-like minima where the vortices 
are strongly localized, the number of local density peaks 
matches the number of vortices N^. The positions of 
these peaks determine what we will call the "vortex lat- 
tice" . 

Much useful information about the spatial structure of 
this vortex lattice may be obtained by plotting directly 
the local density values at the vortex lattice points. An 
excellent, and complementary, alternative to help eluci- 
date the degree of vortex lattice order is to carry out a 
Voronoi construction: we recall that in the Voronoi con- 
struction one determines cells around each lattice point 
by a Wigner-Seitz procedure. In a perfect crystal, all 
the resulting Wigner-Seitz cells are identical, while in 
the general case the cells can have different sizes and 
shapes. The number of sides of the Wigner-Seitz cell 
surrounding a lattice point is identified with the number 
of "nearest neighbors" of that lattice point, and the dif- 
ference between this number and the average (six in our 
case) marks the position of defects (disclinations). We 
will make use of such plots below. 

Since we will be concerned not only about transla- 
tional, but also with orientational order, we will also ex- 
amine the bond-orientational correlation function geir) 
in the vortex lattice defined as: 

geir) = (V'(r)V(O)) (3.2a) 

where the angular brackets denote overall average over 
the vortex lattice and the field ^l^{r) is given by: 

-, "„ 

VXr)^— ^exp[6»0,(r)] (3.2b) 

j=i 

with Oj{r) being the angle that the bond connecting a 
vortex at r to its j-th neighbor makes with a fixed axis, 
and n„ is the number of neighbors of the vortex at r. 

It is also possible, and as we shall see, useful, to define 
what we will call a "translational correlation function" 
in the vortex lattice in a way that is quite analogous to 
the definition of geir). This function, which we denote 
as gci'''): is defined by an equation identical to the right- 
hand side of Eq. H3.2a|l , but with the field being defined 
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FIG. 1: (Color online) The static structure factor S'(k) (see 
Eq. 13. m l for three local minima obtained at T = 18. 4K for 
the same pin configuration. The (red) circles are for a state 
identified, from the overall evidence (see text), as a Bragg 
glass (BrG), the (blue) triangles are for a Bose Glass (BoG) 
state, and the (green) crosses are for an insterstitial liquid 
(IL). The vertical lines are guides to the eye. 



as 

V'G(r) =exp(zG-r) (3.3) 

where G is a reciprocal lattice vector of the triangular 
vortex lattice in the absence of pinning. We will consider 
here only the case where G is a shortest nonzero recip- 
rocal lattice vector and average over the results obtained 
for the six equivalent G's. 

Our identification of the different phases is based on 
analysis of all this information. We show some of the 
results in the next few Figures. These all correspond to 
samples with 64 pins and 4096 vortices. First, in Fig.^ 
we show the structure factor, as defined in Eq. 
The three sets of results shown there are for local free 
energy minima at the same temperature, T — 18. 4K, ob- 
tained with different initial conditions of the three kinds 
described above, and the same pin configuration. Other 
pin configurations yield very similar results (a different 
example is shown in Fig. Ic of Ref. . The green sym- 
bols correspond to the free energy minimum obtained 
by starting with uniform initial conditions. Clearly, the 
structure factor is completely featureless and liquid-like 
in this case - the absolute value of S never exceeds five. 
The red circles are for a local minimum obtained with ini- 
tial conditions corresponding to the best crystalline state, 
as explained above, for this pin configuration. The struc- 
ture factor now exhibits typical ordered behavior, high- 
lighted by six sharp Bragg-like peaks which are empha- 
sized by the added vertical lines in the Figure. Finally, 
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FIG. 2: (Color online). The angularly averaged density cor- 
relation function g{r) plotted at the same temperature as in 
Fig.0as a function of dimensionless distance. The results are 
this time averaged over five different pin configurations. The 
results for the BoG and IL minima are shown in the top panel. 
The (blue) dotted line plot with higher peaks corresponds to 
BoG minima and the (green) light grey plot with fewer and 
lower peaks to IL minima. The bottom panel displays the re- 
sults for the BrG minima. The vertical and horizontal scales 
are the same in the two panels and the color scheme is the 
same as that in Fig.Q 



the blue triangles correspond to a minimum obtained 
by first "quenching" with uniform initial conditions to a 
relatively low temperature (16. 8K) where the liquid-like 
state is found not to be stable, and then slowly warming, 
at 0.2K temperature intervals, back up to 18. 4K. As one 
can see, the >5'(k) for this minimum has an intermediate 
structure, with several relatively well-defined peaks, more 
than six in number, whose heights are considerably lower 
than those of the peaks found for the ordered minimum. 

In Fig. [3 results for the angularly averaged density cor- 
relation function g{r) are displayed. These results are for 
minima obtained at the same temperature, and using the 
same initial condition procedures, as those for the data 
shown in Fig.^ However, the results shown here are av- 



erages over five different pin configurations. The nearly 
flat (green) curve in the top panel corresponds to minima 
obtained from quenches with uniform initial conditions, 
the (red) curve with the well-defined peaks shown in the 
bottom panel corresponds to minima obtained with crys- 
talline initial conditions, and the (blue) curve with the 
intermediate peaks shown in the top panel is for min- 
ima obtained by quenching with uniform initial condi- 
tions to a low temperature and subsequent slow warm- 
ing (very similar results can alternatively be obtained by 
slowly cooling a high-temperature liquid-like minimum 
to a temperature where it is unstable). One can see then 
that the results for g{r) are fully consistent with those 
found in Fig. ^ the minimum obtained from uniform 
initial conditions is fully disordered, while that obtained 
from crystalline initial conditions exhibits a large degree 
of crystalline order. For the the third kind of minima, we 
find some degree of intermediate range order. 

We next examine the structure of the minima in terms 
of the "vortex lattice" found by the procedure explained 
above. Sample results are shown in Fig. O all for the 
same pin configuration, at the temperature T = 18. 2K. 
Results for three minima, obtained from the same proce- 
dures and initial conditions as the three minima in the 
previous two figures, are shown. The positions of the lo- 
cal density peaks and the values of the density at these 
peaks are displayed through a symbol and color coding 
scheme described in the caption of the Figure. Blank 
regions denote areas where no local peak was found us- 
ing the algorithm described above: such areas are mostly 
found only in the most disordered case (uniform initial 
conditions). We see that the results are quite consistent 
with those obtained from the previous correlation func- 
tions: quenching with the appropriate crystalline initial 
conditions leads to a fairly well-ordered (but not per- 
fectly ordered) lattice, while from uniform initial condi- 
tions, one obtains a rather disordered, liquid- like struc- 
ture where the density is nearly uniform and close to 
Pa , except near the pinning centers at each of which one 
vortex is trapped. An intermediate result, on the whole 
more solid than liquid-like, but definitely disordered, is 
obtained through the slow warming or cooling scheme. 

To better gauge the degree of disorder present in each 
case, we construct, in Fig. ^ the corresponding Voronoi 
plots. To display the behavior in a more obvious way, 
the temperature in these plots is lower by 0.2K than 
that of the plots in Fig. |31 As mentioned above, sites 
with number of nearest neighbors (n„) different from six 
represent the locations of disclinations. In the solid-like 
minima (top two panels of Fig. ^J, nearly all the sites with 
n„ 7^ 6 have n„ = 5 or n„ = 7 (the disclinations have unit 
"charge"). Also, six- and seven-coordinates sites, indi- 
cated respectively by (red) triangles and (green) solid cir- 
cles in the plots, always occur in nearest-neighbor pairs. 
Such pairs correspond to dislocations. In the first panel, 
which displays the results for the most ordered state ob- 
tained from crystalline initial configurations, we see that 
the number of dislocations, while nonzero, is quite small 
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FIG. 3: (Color online). Local peak densities at three dif- 
ferent local minima of the free energy. The positions of 
the local density peaks and the values of the density at 
these peaks are displayed according to the following scheme: 
(cyan) solid circles: peaks with pp^^i^/po < 1.5, (red) circles: 

1-5 < Ppeak/Po < 3-5, (blue) squares: 3.5 < /Opeak//^" < 
(green) plus signs: 5.5 < Pp^oi^/po < 7.5, and (purple) tri- 
angles: 7.5 < Ppeak/*"" ^ Bl^nk areas, found mainly in 
the bottom panel, correspond to regions where no local peaks 
are found. The temperature is 18. 2K in all panels, which all 
correspond to the same pin configuration as that for the re- 
sults in Fig0 Pin locations are indicated by (black) asterisks. 
From top to bottom, the results displayed correspond respec- 
tively to ordered (identified as BrG), intermediate (BoG) and 
disordered (IL) minima. 




FIG. 4: (Color online). Voronoi plots for the three cases for 
which density plots are shown in Fig. |3 except that the tem- 
perature is now T — 18. OK. The (black) dots denote lattice 
sites with six neighbors, the (red) triangles denote five-fold 
coordinated sites, and the (green) solid circles, seven-fold co- 
ordinated sites. Rarely occurring four-fold and eight-fold co- 
ordinated sites are indicated by (blue) squares and (purple) 
inverted triangles, respectively. Sites surrounded by black 
circles denote locations of pinning centers. See text for expla- 
nations. 
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and that they form fairly isolated small clusters, each of 
which has zero net Burgers vector. These defect clusters 
occur near the pinning centers and they represent the lo- 
cal disruption of crystalline order due to the pinning. In 
the second panel, which shows the results for the mini- 
mum with intermediate order, the dislocations are orga- 
nized to form well-defined grain boundaries that separate 
crystallites with different orientations. As a result of this 
polycrystalline structure, both translational and bond- 
orientational correlations are short-range. Performing 
the same construction for other pin configurations (re- 
sults for a different pin configuration arc shown in Fig. lb 
of Ref. 13211 ■ one finds that the crystallite arrangement de- 
pends on the pin configuration, as one would expect. The 
grain boundaries tend to lie away from pinning centers, 
which makes sense physically. Finally, in the last panel, 
which shows the results for the minimum obtained from 
uniform initial conditions, the high degree of disorder is 
evident. The Voronoi construction shown in this panel 
is not very meaningful because our method of obtaining 
the vortex positions is not reliable for liquid-like minima. 
The identification of the position of a vortex with that of 
a local peak of the density field is justified only when the 
peak is sharp. This is the case in the solid-like minima, 
but not so in the liquid-like minima with low peak den- 
sities. We have shown here the results for the liquid-like 
minimum only for the purpose of illustrating the differ- 
ence between the structure of this minimum and those of 
the more ordered ones. For this minimum, we find a very 
large number of defects, and it is difficult to determine 
whether disclinations of opposite "charge" always pair up 
to form dislocations. Inspection of the defect distribution 
suggests the presence of free disclinations, but not con- 
clusively. The total number of local peaks of the density 
field in liquid-like minima turns out to be substantially 
smaller than the expected number of vortices. This prob- 
lem is not present at the solid-like minima where nearly 
all the vortices are strongly localized. 

Next, we show in Fig. \E\ examples of the bond- 
orientational and "translational" correlation functions, 
geir) and gci'r) as defined in Eqs. (|3.2|l and H3.3|) respec- 
tively. All results shown in this figure are averages over 
five pin configurations at temperature T = 17. 6K. At this 
temperature the liquid is unstable, and in any case all liq- 
uid correlation functions are featureless, so that case is 
not shown. The (purple) triangles correspond to ordered 
minima obtained from quenching with crystalline initial 
conditions to T = 16. 8K and subsequent slow warming, 
and represent the bond-orientational correlation, gQ{r). 
This function appears to saturate at a fairly large value 
as r increases, indicating the presence of long-range bond- 
orientational order. The (red) open circles display g(,{r) 
for the polycrystalline state with intermediate order, ob- 
tained through the procedures described above. This 
function decays exponentially with distance, as shown 
by the exponential fit (dashed line), over a distance of 
few tens in units of ao. Finally, the (black) open squares 
show the results for the "translational" correlation func- 
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FIG. 5: (Color online). Bond-orientational correlation func- 
tion gair) as defined in Eq. 113. 2|l . and "translational" correla- 
tion function, gG{r), see Eq. 13.311 . These are results averaged 
over five pin configurations, at T = 17. 6K. The (purple) tri- 
angles and (red) circles denote gsir) for ordered (BrG) and 
polycrystalline (BoG) minima (liquid-like minima are not sta- 
ble at this lower T) respectively. The black squares are gcif) 
for the BrG minima. The dashed line is an exponential fit to 
g(i{r) for the BoG minima. The dotted line is a similar ex- 
ponential fit to ga{r) for the BrG minima at small distances, 
showing that the decay is slower than exponential for large 
values of r. The solid lines connecting symbols are guides to 
the eye. 



tion go (r) for the ordered minima. The dotted line is an 
exponential fit to this translational correlation function 
for r/flo < 10, showing that the decay at longer distances 
is clearly slower than exponential (possibly a power-law). 

^From the overall examination of data such as those 
shown in these figures, we can reach the following con- 
clusions: quenching with uniform initial conditions to 
temperatures above or slightly below leads to very 
disordered states, with very short range correlations and 
no structure, except for vortices at and near the pin- 
ning centers. These minima become unstable upon slow 
cooling: one then obtains minima of the third kind, dis- 
cussed below. They clearly must be identified with the 
liquid phase, specifically the interstitial liquid (IL) phase 
discussed in the Introduction. The second kind of min- 
ima are obtained by quenching with crystalline initial 
conditions to any temperature below, or slightly above, 
T^. These minima can then be cooled, without losing 
their character, but they melt into the IL upon suffi- 
cient warming. The states thus obtained are obviously 
nearly crystalline, with the reservation that the order is 
not truly long range, but has a slow decay. This is most 
evident in the results for gair) seen in Fig. [3 The bond- 
orientational order is nearly perfect. Defects are limited 
to isolated clusters. This state must therefore be identi- 
fied as a Bragg glass (BrG). Finally, we have the states 
obtained as described above, by slow cooling of the IL 
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state, which can also be obtained by quenching with uni- 
form initial conditions to a temperature below the IL 
state stability limit. Either procedure yields, for any pin 
configuration, states that at a given T differ only slightly 
in free energy or density configuration. The states do dif- 
fer somewhat more, however, for different pin configura- 
tions. The Voronoi construction conclusively shows (see 
Fig. ^ that these states are polycrystalline. The grain 
boundaries are formed by dislocation chains. For any 
pin configuration, these grain boundaries survive thermal 
cycling across T^, but if sufficiently warmed up, these 
minima melt, the melting beginning locally at the grain 
boundaries (see Sec lIII BI below'). The bond-orientational 
correlation function for these minima has an intermediate 
range. This, and also the observation that these states 
are not unique (as mentioned above, slightly different 
states are reached depending on the cooling or quench- 
ing protocol), clearly indicate a glassy state. We there- 
fore identify it with the Bose glass (BoG) which has been 
experimentalljiiSiiS and numericalljiS^ shown to be poly- 
crystalline. 



B. Properties of the condensed phases 

In this section, we describe in detail how the proper- 
ties of the low-temperature "solid-like" condensed phases 
represented by the BrG and BoG types of minima de- 
fined in the preceding section vary as the temperature is 
changed. The high-temperature liquid phase represented 
by the IL-type minima is not very interesting from the 
point of view of its temperature dependence: in the tem- 
perature range we have considered, this phase exhibits 
a liquid-like (nearly uniform) density distribution except 
in the immediate vicinity of the columnar pins, at each of 
which a vortex is trapped. These trapped vortices would 
eventually get delocalized at higher temperatures. This 
happens beyond the upper limit of the temperature range 
(16.8K < T < 18. 8K) we have considered here. 

To study how the properties of a local minimum of 
the free energy change as the temperature is varied, we 
have "followed" minima as the temperature changes: for 
example, starting with a minimum obtained by quench- 
ing to the lowest temperature of 16. 8K, we "follow" that 
minimum to higher temperatures by increasing the tem- 
perature in small steps (usually taken to be 0.2K) and 
finding a new minimum at the higher temperature by 
running the minimization routine with the configuration 
at the minimum obtained at the previous temperature as 
the initial state. This procedure leads to a new minimum 
of the same type as long as the minimum remains locally 
stable - if the temperature is increased to values substan- 
tially higher than T^, where the BrG and BoG minima 
become unstable, this procedure leads then to the IL min- 
ima. We have also carried out "cooling" runs where the 
temperature was decreased in small steps, starting from 
a minimum obtained at a relatively high temperature. 
When the IL minima become unstable, BoG states are 
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FIG. 6: (Color online). Temperature dependence of the order 
parameter m as defined in Eq. (13.411 for three kinds of minima, 
averaged over five samples, with Nv = 4096 and N-p — 64. 
Results for the BrG, BoG and IL minima are shown by (red) 
plus signs, (green) crosses and (blue) asterisks, respectively. 

obtained. As long as the minimum under consideration 
did not become unstable in the range of temperatures 
considered in the runs, we did not find any substantial 
difference between the results obtained in the "heating" 
and "cooling" runs 

To characterize the density distribution at a minimum, 
we looked at "global" quantities such as the structure fac- 
tor S'(k), an example of which is shown above in Fig.^ 
The T-dependent information contained in 5'(k) is, how- 
ever, difficult to display in plots. We therefore introduce 
a closely related "order parameter" m defined as 

m = ^{S^ax/N^) (3.4) 

where Smax is the largest value of the structure factor 
S'(k) averaged over the six k- vectors related by lattice 
symmetry to each k. By definition, the order parameter 
m is equal to unity in a state with perfect crystalline or- 
der (triangular lattice with (5- function peaks). Since no 
long-range crystalline order is expected or found in either 
one of the BrG and BoG phases, the value of m should go 
to zero at all temperatures in the thermodynamic limit. 
However, the infinite size limit is reached very slowly in 
the glassy phases, and for the finite-size systems consid- 
ered here, this quantity provides a convenient measure of 
the degree of local order present at the minima. 

We show in Fig. the results for the order parame- 
ter m, averaged over five different pin configurations, for 
samples with 4096 vortices and 64 pins. The results for 
the BrG, BoG and IL minima are shown by (red) plus 
signs, (green) crosses and (blue) asterisks, respectively. 
Only two data points are shown for the IL minimum be- 
cause it becomes unstable at lower temperatures. For 
all three kinds of minima, the value of m decreases with 
increasing T, as expected. The rate of change is largest 
for the BrG minima, and smallest for the IL minima. 
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FIG. 7: (Color online). Peak density plots for the BoG min- 
imum of a sample with — 4096 and Np = 64 (the sample 
is the same as the one for which results are shown in Figs. |21 
and^ at two temperatures: T = 17.4K (top panel) and T — 
17.8K (bottom panel). The symbols and color scheme used 
here are the same as those in Fig. |3 

However, the differences among the three kinds of min- 
ima in the degree of local order remain quite clear at 
all the temperatures considered. We have also obtained 
BrG-type minima for several (~ 10) samples with 1024 
vortices and 16 pins. The values of m obtained from an 
average over these smaller samples are found to be only 
about 6% larger than those obtained for the 4096- vortex 
samples at the corresponding temperature. This result 
implies that the translational correlation function at the 
BrG minima falls off very slowly with distance. This 
is consistent with our interpretation of these minima as 
representing the BrG phase. 

We also examined the temperature dependence of the 
detailed density distribution at the minima (obtained 
from "peak density" plots similar to those shown in 
Fig. I^l and their defect structure obtained from Voronoi 
plots such as those in Fig. 01 Examination of the peak- 
density plots reveals more information about how the 



density distribution at the minima changes with temper- 
ature. In Fig. [7| we have shown peak-density plots at 
two temperatures (17. 4K in the top panel and 17. 8K in 
the bottom panel) for the BoG minimum of the same 
sample for which a similar plot at 18. 2K is shown in the 
middle panel of Fig. O The symbols and color scheme 
used in these plots are the same as those in Fig. O It is 
clear from these plots that the typical values of the local 
peak densities decrease (the vortices become less local- 
ized) as the temperature increases. The largest values 
of the local peak density occur near the pinning centers 
and the smallest values appear near the grain bound- 
aries that separate different crystalline regions. This is 
physically reasonable - the disruption of local crystalline 
order near the grain boundaries should make the vor- 
tices in such regions more delocalized. At the relatively 
low temperature of 17. 4K, there are very few local den- 
sity peaks where the peak density is lower than 3.5po 
(such peaks are indicated by (red) circles in the plot). 
As the temperature is increased to 17. 8K, the number 
of such peaks increases, and this trend continues as T is 
increased further, as can be seen in the middle panel of 
Fig. O where the results for T = 18.2K are shown. It is 
also clear from these plots that the spatial regions where 
the low peak densities occur remain roughly unchanged 
as the temperature increases - as noted above, these re- 
gions are strongly correlated with the locations of the 
grain boundaries. These observations indicate that the 
vortices in the interior of the crystalline grains remain 
in a "solid" state as the temperature is increased to a 
value close to T^, while those in the neighborhood of the 
grain boundaries begin to "melt" (get delocalized) at a 
lower temperature. This "inhomogeneity" of the melting 
process will be discussed in more detain in section UlI Dl 

We also studied the temperature dependence of the 
number of topological defects present at the minima. 
This number increases slowly with increasing T, and then 
exhibits a sudden jump as liquid-like regions (character- 
ized by low values of the local peak densities) appear 
near the grain boundaries at temperatures close to . 
As noted in Sec. lIII Al above. our method of obtaining the 
defect structure using the Voronoi construction becomes 
somewhat less reliable when liquid-like regions begin to 
appear. Due to this difficulty in obtaining reliable re- 
sults for the defect structure and statistics at relatively 
high temperatures, we have not carried out a quantita- 
tive analysis of the dependence of these quantities on the 
temperature. 

The gradual decrease in the values of the local peak 
densities with increasing temperature is also found in 
the BrG minima, but on the average the vortices remain 
more strongly localized at the BrG minima than at the 
corresponding BoG ones. The strong correlation between 
low values of the local peak density and the location of 
topological defects is found at the BrG minima too. The 
main difference between BrG and BoG minima is that 
the defects are more randomly distributed at the BrG 
minima: they do not line up along grain boundaries as 
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they do at the BoG minima. For this reason, the regions 
of low peak density appear at fairly random locations at 
the BrG minima. The total number of topological de- 
fects at the BrG minima remains smaller than that at 
corresponding BoG minima at all temperatures. As the 
temperature approaches T^, liquid-like regions appear 
in parts of the sample where the defect density is large. 
This is illustrated in the plots in the top panels of Figs.|3| 
and|31 

Thus, our investigation of the temperature dependence 
of the structure of the BrG and BoG minima confirms 
that these minima retain their topological structure as T 
is increased toward and even beyond (provided they 
remain stable). The degree of localization of the vortices 
decreases with increasing T and liquid-like regions with 
nearly delocalized vortices begin to appear at tempera- 
tures close to T^. However, even at such temperatures, 
the BrG and BoG type minima are clearly distinguish- 
able from one another and from the IL-type minima, as 
illustrated in Figs |31 E and El 



C. Two-step melting transition 

As discussed in the preceding subsections, we find three 
kinds of coexisting, locally stable minima of the free en- 
ergy at temperatures close to T^. Two kinds (BrG and 
BoG) of local minima are found at temperatures substan- 
tially lower than T,^, and only the IL minimum is stable 
if T is much higher than T^. In our mean-field descrip- 
tion, the thermodynamically stable phase at a particu- 
lar temperature corresponds to the minimum with the 
lowest free energy. Therefore, crossings of the free ener- 
gies of different kinds of minima correspond to first-order 
phase transitions in our description. The IL minima are 
expected to have the lowest free energy at high tempera- 
tures and the solid-like BrG or BoG minima should rep- 
resent the globally stable phase at low temperatures. To 
determine how this "freezing" transition takes place, it is 
necessary to examine the dependence of the free energies 
of these different kinds of minima on the temperature T. 

In Fig.|Sl we have shown the results for the free energies 
of different minima of a sample with Ny — 4096 and 
Np = 64. This sample is the same as the one for which 
results are shown in Figs. ^ El El and[7| A similar plot, 
showing the same general behavior for a different sample, 
may be found in Ref. [33 . The BrG minimum clearly 
has the lowest free energy at low temperatures (while 
we have shown data for T > 17.8K, we have checked 
that this remains true at lower temperatures), while the 
IL minimum is the one with the lowest free energy at 
high temperatures. The interesting feature that emerges 
from data as that shown in this Figure is that the BoG 
minimum has the lowest free energy in an intermediate 
temperature range of small width (between 18. OK and 
18. 3K for this sample), indicating that the melting of the 
low-temperature solid phase to the high-temperature IL 
phase occurs in two distinct steps: the BrG phase that is 
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FIG. 8: (Color online). Temperature dependence of the di- 
mensionless free energies of different local minima of a sam- 
ple with — 4096 and Np — 64 (the sample is the same as 
the one for which results are shown in Figs. El El ''■iid [7J. 
The data for the BrG, BoG and IL minima are shown by the 
(green) crosses and dashed line, (blue) asterisks and dotted 
line and the (red) plus signs and solid line, respectively. The 
lines are guides to the eye. 



the thermodynamically stable one at low temperatures 
undergoes a first-order transition into a BoG phase as 
the temperature is increased, and then this BoG phase 
melts into the IL phase via a second first-order transition 
at a slightly higher temperature. This two-step melting 
behavior is one of our main results. 

The same qualitative behavior, indicating the occur- 
rence of two separate first-order transitions, is found in 
all the 4096-vortex samples we have studied at pin con- 
centration c = 1/64 {Np = 64). The average value of 
the temperature interval in which the BoG phase has the 
lowest free energy is about 0.42K. This width exhibits 
fairly large sample-to-sample variations, ranging between 
about O.IK in one sample to a maximum of 1.2K. As the 
system is cooled from the high-temperature liquid phase, 
it undergoes a first-order transition into a polycrystalline 
BoG phase. Such a transition has been observedi^ in ex- 
periments on BSCCO samples with a small concentration 
of columnar pins. The value of the upper (BoG to IL) 
transition lies between 18. 2K and 18. 3K in all the 4096- 
pin samples we have studied. These values are quite 
close to the first-order melting temperature of the 
same system in the absence of pins'^^ . Our results, thus, 
are consistent with the experimental observationiii^ of a 
weak dependence of the freezing temperature of the vor- 
tex liquid on the pin concentration c for small values of c. 
In addition, our work predicts a second first-order tran- 
sition to a more ordered BrG phase at a slightly lower 
temperature. 

For smaller samples at the same c = 1/64 concentra- 
tion {Ny = 1024 with Np = 16) we did not always find the 
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FIG. 9: (Color online). Results for a sample with Nv = 1024 
and Np = 16. Voronoi plots for the BrG and BoG minima at 
18. OK are shown in the top and middle panels, respectively. 
The symbols and color scheme used in these plots are the same 
as those in Fig. |1] The bottom panel shows the temperature 
dependence of the free energies of the BrG (blue line and 
squares), BoG (red line and circles) and IL (black line and 
triangles) minima, respectively. The lines are guides to the 
eye. 
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FIG. 10: (Color online). The pinning free energy (see text) 
averaged over five samples with = 4096 and Np — 64. 
Results for the BrG and BoG minima are shown by (red) 
triangles and (black) circles, respectively. The solid lines are 
straight lines joining the points. 



BoG minima: the minimization procedure that led to the 
BoG minima in the samples with — 4096 often con- 
verged to minima of the BrG type in the smaUer samples. 
This is because the typical size of the crystalline grains at 
the BoG minima at this pin concentration is of the order 
of the sample size for 1024-vortex samples. This can be 
seen by comparing the middle panel of Fig. ^ where we 
have shown the Voronoi plot for the BoG minimum ob- 
tained for a A^^ = 1024 sample, with the Voronoi plot in 
the middle panel of Fig. 01 which is for the same T and c, 
but a larger sample with 4096 vortices. One can plainly 
see that the domain size is comparable to the system 
size of the smaller samples. However, for the 1024-vortex 
samples where we found BoG-type minima, the behav- 
ior of the free energies was found to be very similar to 
that shown in Fig. |S| An example of such behavior is 
shown in the bottom panel of Fig, |51 The two-step melt- 
ing transition found in the larger systems is found here 
also, indicating that this is the generic behavior. We 
have also carried out preliminary studies of the phase 
diagram for a slightly higher value of the pin concentra- 
tion, c = 1/32. We find the same qualitative behavior, 
namely the occurrence of two-step melting, with a larger 
difference between the two transition temperatures. We 
therefore conclude that the two-step melting transition 
we have found is a characteristic feature of our model 
of layered superconductors with a small concentration of 
columnar pins. 

Since a first-order transition is not characterized by 
power-law behaviors of thermodynamic quantities with 
universal exponents, the usual finite-size scaling analysis 
of numerical data for continuous transitions does not ap- 
ply to our work. However, we have examined the depen- 
dence of the magnitudes of various discontinuous changes 
at the transitions on sample size. From the results for 
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the free energy as a function of the temperature, we have 
calculated the entropy jump per vortex at the first-order 
transitions. The entropy change at the BrG-BoG (lower 
temperature) transition is found to be ~ O.lfcs per vortex 
{ks is the Boltzmann constant), and the corresponding 
change at the BoG-IL (higher temperature) transition 
has a slightly higher value, ~ 0.15^^ per vortex. The 
sum of these two entropy jumps is slightly smaller than 
the entropy change (~ 0.29fc_B per vortex) at the single 
first-order melting transition found^l in the same system 
in the absence of pinning. All these results are physi- 
cally reasonable. The difference between the degrees of 
order in the BrG and BoG minima is smaller than that 
between the BoG and IL minima, suggesting that the 
entropy jump at the lower transition should be smaller 
than that at the upper one. Also, the low-temperature 
BoG phase of the system in the presence of pinning is 
less ordered than the crystalline phase of the pure sys- 
tem, and the IL phase is slightly more ordered than the 
vortex liquid in the system without pinning (this is due 
to the local order of the vortices near the pinning cen- 
ters which persists in the IL phase). Therefore, the net 
change of entropy in going from the BoG phase to the 
IL phase should be smaller than the entropy jump at the 
melting transition in the pure system. A comparison of 
our results for two sample sizes {N^ — 4096 and 1024) 
does not show any appreciable dependence of the entropy 
changes on the sample size. We, therefore, conclude that 
there is no indication of any significant change in our 
phase diagram as the sample size is increased. 

Since our mean-field treatment ignores the effects of 
fluctuations, it is important to address the question of 
whether one or both of the transitions found would be- 
come continuous if fluctuations were included. However, 
there are very few examples of such fluctuation-driven 
continuous transitions in three dimensions. In our cal- 
culations, the effects of the electromagnetic interaction 
among vortices on different layers are included. So, we 
expect our results regarding the nature of the transi- 
tion to remain valid if fluctuations were included. This 
conclusion is supported by the result that our density 
functional theory provides a quantitatively a correct ac- 
(,Qyj^-^3ij33j34 -j-j^g first-order melting of the vortex lattice 
in the pure system. Also, the transition in the presence 
of a small concentration of columnar defects is known to 
be first-order both experimentalljjiiiSiiS and from simu- 
lations^. 

The appearance of a sliver of the BoG phase in our 
phase diagram may be qualitatively understood as aris- 
ing from a competition between the elastic and pinning 
parts of the free energy. The elastic (free) energy plays 
a dominant role in the BrG minima: the effect of the 
randomly located pinning centers is accommodated in 
this structure by small displacements of the vortices from 
their ideal lattice positions towards the nearest pinning 
center. Where the occupation of a pinning center by a 
vortex would require a large displacement of the vortex, 
the pinning center is not fully occupied. For this reason. 



the number of pinning centers occupied by vortices at 
a BrG minimum is always slightly lower than the total 
number of pinning centers. At the corresponding BoG 
minima, on the other hand, the pinning centers are al- 
ways occupied (except in very rare cases where two pin- 
ning centers are located very close to each other). Since 
the BoG minima are obtained from liquid-like initial con- 
ditions, the starting point has fully occupied pinning cen- 
ters, and vortices localized in small crystalline patches in 
their immediate vicinity (see, for example, the bottom 
panel of Fig. O which shows such crystalline patches sur- 
rounding the pinning centers). The crystalline orienta- 
tion of vortices around a vortex-binding pin site depends 
on the local arrangement of the pinsi2& This orientation 
is, in general, different in different regions of the sample. 
These crystalline patches grow as the temperature is re- 
duced, and eventually meet one another at grain bound- 
aries to form a polycrystalline BoG minimum. It is clear 
that the pinning centers are better accommodated in this 
structure than in the corresponding BrG structure. On 
the other hand, the creation of grain boundaries costs 
elastic energy. Thus, the BoG minima are expected to 
have lower pinning energy but higher elastic energy than 
the BrG ones for the same pin configuration. The elastic 
energy dominates over the pinning energy at low tem- 
peratures in our low pin-concentration samples. This is 
why the BrG phase is globally stable at low T. Both 
these components of the free energy decrease in magni- 
tude as the temperature is increased. The softening of 
the lattice near the melting transition causes the elastic 
energy to decrease faster than the pinning energy. This 
makes the pinning component of the free energy more 
important than the elastic part near T^, thereby mak- 
ing the total free energy of the BoG minima (which, as 
discussed above, have lower pinning energy) lower than 
that of the BoG minima. The BoG to IL transition at 
a slightly higher temperature is driven by the usual en- 
tropic mechanism. 

We can substantiate this qualitative explanation by 
evaluating the pinning component of the free energy, 
as defined in Eq. H2.4|) . for the BrG and BoG minima. 
The results of our calculation, averaged over five sam- 
ples (each with 4096 vortices and 64 pins), are shown in 
Fig. 1101 The pinning energy is negative for both kinds of 
minima, as expected. The absolute value of the pinning 
energy of the BoG minima shows a smooth (nearly linear) 
decrease as the temperature is increased. This is because 
both the value of the parameter F (see Eq. H2.6|l ) that de- 
termines the depth of the pining potential, and the height 
of the local density peak at a pinning center decrease with 
increasing T. The plot for the BrG minima shows a sim- 
ilar but much slower behavior, still monotonic or nearly 
so: it is not clear whether the very shallow minimum 
near T = 17. 6K is significant. In some of the samples, 
the total number of occupied pining centers increases by 
a small amount near this temperature, thereby reducing 
the value of the total pinning energy. This may reflect a 
better accommodation of the pinning centers by the lat- 
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tice, which softens as the temperature is increased. This 
plot clearly demonstrates that the pinning centers are 
better accommodated at the BoG minima. The differ- 
ence in the pinning energies of the two kinds of minima 
decreases with increasing T, but the difference in the elas- 
tic component of the free energy decreases faster (this is 
clear from our results for the total free energy) due to a 
softening of the elastic constants with increasing temper- 
ature. The overall effect is a crossing of the free energies 
of the two kinds of minima near the melting transition. 

Quantitatively, however, the values of the two parts 
of the free energy and their dependence on T are deter- 
mined by the material parameters of the superconductor 
and the properties of the pinning centers. Since several 
experimental studiesiiiiSiiS of the effects of irradiation- 
induced columnar pinning on the mixed state of BSCCO 
exist in the literature, we have used parameter values 
appropriate for this system. For other layered super- 
conductors, the BoG phase might occur over a wider or 
narrower (even vanishing) range. 



D. Inhomogeneous melting of the vortex solid 



As mentioned in section Fill Bl the density distribution 
in both BrG and BoG minima is very inhomogeneous: 
at temperatures close to Tj^ there are liquid-like regions, 
characterized by low values of the local peak densities 
in parts of the sample. In this subsection, we show that 
this inhomogeneity leads to a spatial variation of a "local 
melting temperature" , defined below. As mentioned in 
the Introduction, such spatial variations of a local melt- 
ing temperature have been deduced from measurements 
of the local magnetization in BSCCO samples with ran- 
dom point^^ and columnar pinningiS. 

Since we are dealing with time averaged densities, it is 
clear that the value of pi at a local peak of the density 
field provides a measure of the degree of localization of 
the vortex whose average position corresponds to the lo- 
cation of the density peak. A high (low) value of the local 
peak density implies strong (weak) localization. Smaller 
values of the local peak density imply mobile, liquid-like 
behavior. The value of the local density is, of course, 
equal to po everywhere in the liquid state in the absence 
of pinning. The presence of pinning centers causes the 
local density to vary in space, but in the liquid state this 
variation does not lead to local peaks higher than about 
Spo (excluding the vortices localized at the pinning cen- 
ters) . A similar result was also obtained in previous stud- 
ies^ of vortices in the presence of pinning, where it was 
found that values of the peak density lower than about 
Spo correspond to the liquid state. We, therefore, take 
the value 3po of the local peak density as separating solid- 
and liquid-like behaviors. 

Using this criterion, we can determine whether a small 
region of the sample at a given minimum is in a locally 
"solid" or "liquid" state. To do this, we define a quantity 
as the average of the local peak densities in a small 
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FIG. 11: (Color online). Determination of the local melting 
temperatures in different regions of a sample. The top panel, 
at T = 18. OK and for the BoG minimum, shows the three 
different regions (A, B and C) considered in the calculation. 
The positions of local density peaks in these regions are in- 
dicated by (red) plus signs (region A), (blue) crosses (region 
B) and (green) circles (region C). The black dots represent 
the positions of the local density peaks in the other regions, 
and the (black) circles show the locations of the pinning cen- 
ters. The bottom panel shows the dependence of the average 
local peak density p^v (normalized by the liquid density po) 
calculated for these regions on the temperature T. Data for 
regions A, B and C are respectively shown by (red) triangles, 
(blue) squares, and (green) circles. The solid lines are guides 
to the eye. 



region containing ^ 100 vortices (see Fig. Illf) . The very 
high local density peaks representing vortices trapped at 
pinning centers are not included in this average. Values 
of pP„ substantially larger than 3po indicate solid-like be- 
havior in the region under consideration, while substan- 
tially lower values of p^^, suggest a locally melted region. 
To determine how the local melting temperature varies 
from one part of the sample to another, we have studied 
the temperature dependence of p^^ for different regions 



of the sample. Typical results are shown in Fig. for 
the pin configuration for which we have earlier shown de- 
tailed results in several Figures. In the top panel, where 
the positions of the local density peaks in the BoG min- 
imum at T = 18. OK and the locations of the pinning 
centers are shown, we have indicated three regions for 
which the average local peak density pg^ was calculated. 
The local density peaks in these regions are indicated by 
(red) plus signs (region A), (blue) crosses (region B) and 
(green) circles (region C), while the black dots represent 
the positions of the local density peaks in the other re- 
gions. The bottom panel of Fig.llllshows the dependence 
of the p^yS calculated for these three regions on the tem- 
perature T . At each value of T, the minimum with the 
lowest free energy at that temperature was used in the 
calculation of p^^. 

It is clear from the plot and from similar ones for other 
pin configurations (one of which was shown in Ref. 1321) 
that the "local melting temperature" , operationally de- 
fined as the temperature at which drops from values 
well above 3po to values clearly below, varies from region 
to region. The range over which the local melting tem- 
perature varies is comparable to that found in the exper- 
imentsiSiS^. The values of the local melting temperature 
are strongly correlated with the local pin structure and 
the resulting defect structure of the vortex solid. Region 
B of the sample does not have any pinning center and is 
located across a grain boundary of the BoG minimum. 
As noted in section Fill Bl regions near grain boundaries 
appear to melt at temperatures that are slightly lower 
than the global transition temperature determined by the 
crossing of free energies. This is why the local melting 
temperature for region B is the lowest. This tempera- 
ture corresponds to the transition from the BrG to the 
BoG state. In region A the vortices form a nearly perfect 
crystalline arrangement with no topological defects. The 
local melting temperature measured in this region is, as 
a consequence, close to that of the pure system. Region 
C contains a cluster of pinning centers, and the relatively 
large value of for this region at temperatures higher 
than the BoG to IL transition temperature reflects the 
local solid-like structure of vortices situated near pinning 
centers. Thus, the spatial variation of the local transi- 
tion temperature is closely correlated with the "pinning 
landscape" of the sample. 

We have found a second way of correlating the melt- 
ing transition with the local density structure of the free 
energy minima. This second way is based on a previ- 
ous study^^ of the melting of the vortex lattice in the 
presence of periodic pinning, where it was found that the 
melting transition corresponds to the onset of percolation 
of liquid-like regions defined using the peak-density crite- 
rion mentioned above. We focus for this purpose on the 
upper transition between the BoG and IL phases. We 
classify vortices represented by local peaks of the den- 
sity as solid-like or liquid-like, depending on whether the 
density at the local peak is higher or lower than 3po- 
We then check whether the regions containing liquid-like 
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FIG. 12: (Color online). Percolation and melting (see text). 
The (red) plus signs denote the locations of liquid-like local 
density peaks (peak height smaller than 3po), and the (blue) 
circles indicate the solid-like ones (peak height greater than 
3po). The top panel is for the BoG minimum of a 4096- vortex 
sample at 18. 2K and the bottom panel is for the IL minimum 
of the same sample at 18. 4K. The thermodynamic transition 
between the BoG and IL phases takes place between these 
two temperatures. 



vortices percolate across the sample. We find, as in the 
earlier study, that the melting transition coincides with 
the occurrence of percolation of the liquid-like regions. 
Typical results are shown in Fig. 1121 where the top panel 
shown the locations of solid- and liquid-like local density 
peaks at a BoG minimum obtained at 18. 2K (this mini- 
mum has the lowest free energy at this temperature) , and 
the bottom panel shows a similar plot for the IL mini- 
mum of the same sample at 18. 4K (the IL minimum is 
the one with the lowest free energy at this temperature, 
so that the melting transition in this sample occurs be- 
tween 18. 2K and 18. 4K). The liquid-like regions do not 
percolate at 18. 2K, whereas they do at the higher tem- 
perature of 18.4K. Thus, the melting transition in this 
sample also corresponds to the occurrence of percolation 
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for the liquid-like regions. Similar behavior was found 
in the other samples we have studied. This observation 
provides a convenient way of approximately locating the 
transition point. This may be useful in other situations 
where the method of locating the transition temperature 
from a crossing of free energies may not be easily imple- 
mentable (e.g. at higher pin concentrations where the 
melting transition is expected to become continuous). 



IV. CONCLUSIONS 

We have presented here the results of a detailed inves- 
tigation of the structural and thermodynamic properties 
of a system of vortices in a highly anisotropic layered 
superconductor with a small concentration of randomly 
placed columnar pinning centers. Both the external mag- 
netic field and the columnar pins are assumed to be per- 
pendicular to the superconducting layers. Our method, 
based on numerical minimization of the appropriate free- 
energy functional, allows us to obtain very detailed in- 
formation about the density distribution in the different 
free energy minima (which allows us to reliably identify 
the phases corresponding to these minima), and to map 
out the phase diagram of the system. 

There are several salient results of our study. The 
first one is the occurrence of a topologically ordered BrG 
phase at low temperatures. While the occurrence of 
a low-temperature BrG phase in superconductors with 
a low concentration of random point pinning centers is 
well-established now, relatively little is known about the 
existence of such a phase in superconductors with random 
columnar pins. Our results are consistent with those of a 
recent numerical studjJ^S, of a similar system. It is clear 
from our work that the BrG minima represent a phase 
distinct from the polycrystalline BoG phase also found in 
our study. We can not conclusively rule out the possibil- 
ity that free dislocations would appear at the nearly crys- 
talline minima at length scales much longer than those 
considered in our numerical study. If this should happen, 
then the "hexatic glass" phase suggested in some earlier 
theoretical studies"^^ would become a possible candidate 
for describing the BrG minima found in our study. This 
phase, however, would be distinct from the BoG. 

Our second important result is the occurrence of a two- 
step melting transition: we find that the low-temperature 
BrG phase transforms into a polycrystalline BoG phase 
as the temperature is increased, and this BoG phase then 
melts into the high-temperature IL at a slightly higher 
temperature. This conclusion about the occurrence of 
two distinct transitions would remain valid even if the 
true nature of the low-temperature phase turns out to 
be slightly different from a Bragg glass: our work shows 
that the BrG and BoG minima are quite distinct from 
each other. The possibility of occurrence of a two-step 
melting transition of the vortex lattice in systems with 
random point pinning has been suggested earlier—. Our 
work provides support to this suggestion. 



To our knowledge, the second (lower T) transition be- 
tween the BoG and BrG phases has not been observed 
in experiments on layered superconductors with a small 
concentration of random columnar pins. This may be due 
to strong metastability: the BoG minimum into which 
the IL is expected to freeze as the temperature is de- 
creased remains locally stable at temperatures lower than 
that at which its free energy crosses that of the BrG mini- 
mum, suggesting that it would be difficult to see in exper- 
iments the transition to the globally stable BrG phase. 
The situation here may be similar to that found in a 
recent experimental study^ of a low- Tc superconductor 
with weak point pinning which is expected to exhibit a 
BrG phase at low T. It is, however, found in the ex- 
periment that the vortex solid obtained by cooling the 
sample in the presence of the external magnetic field has 
a polycrystalline structure, indicating that the metasta- 
bility of this disordered state prevents the system from 
reaching the more ordered (BrG) equilibrium state at low 
temperatures. In our numerical work, the BrG minima 
were obtained by performing the free-energy minimiza- 
tion from an initial configuration with crystalline order. 
It is not clear how a similar procedure can be adopted in 
experimental studies. 

A low-temperature BrG phase has been observed in 
a recent simulationSi of a model of HTSC with a small 
concentration of columnar pinning centers. However, this 
simulation finds a single first-order transition between the 
BrG and IL phases at low pin concentrations: the small 
intermediate region of BoG phase found in our study is 
not observed. This is probably due to the smallness of the 
system sizes (~ 100 vortex lines) used in the simulation. 
As discussed in Sec. lIII Cl above. the polycrystalline BoG 
minima are found only if the sample size is larger than 
the typical size of the crystalline domains: only the BrG 
minimum is found at low temperatures if this condition is 
not satisfied. Since the crystalline domains become large 
at small pin concentrations, it is quite likely that a simu- 
lation with small system sizes and low pin concentration 
would not see the BoG phase. Another possibility is that 
a narrow "two-phase" region found in Ref. near the 
BrG melting transition actually corresponds to the inter- 
mediate BoG phase found in our study. A third possibil- 
ity is that the intermediate BoG phase, which arises due 
to a competition between the elastic and pinning compo- 
nents of the free energy (see Sec. 1111 01 where it is shown 
that the pinning energy favors the BoG phase), does not 
appear in the simulation of Ref. |^ because the pins are 
assumed to be "weak" in that work. 

Our study also illustrates the spatial inhomogeneity of 
the melting process in the presence of disorder: we have 
found that a "local" melting temperature defined using 
a criterion based on the degree of localization of the vor- 
tices shows considerable spatial variation. This spatial 
variation is strongly correlated with the "pinning land- 
scape" associated with the random spatial location of the 
pinning centers. Our results about the spatial inhomo- 
geneity of the local melting temperature are consistent 
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with those of recent experiments^^^ on superconductors 
with pinning disorder. 

These results estabhsh the usefulness of our numerical 
method in dealing with the problem of vortex matter in 
the presence of random pinning. It would be interesting 
to examine the extent to which our results depend on 
the values of the parameters appearing in the free en- 
ergy functional. It would also be clearly useful to carry 
out similar studies of other related problems, such as the 
complete phase diagram of systems with random colum- 



nar pins in the T — c plane, and the phase behavior of 
systems with random point pinning. Some of these in- 
vestigations are currently in progress. 
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